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1.  INTRODUCTION 


Let  X  >  (X,,...,X  )'  be  a  nonsingular  multivariate  normal  random  vector  which  is 

1  n 


standardized  so  that  E(X)  *  0  and  var(X)  «  R,  where  R  *  (/>..)  with  p  *  1  for  all  i,  and 

9  II 


let  A«Rn  be  e  permutation-symmetric  region.  In  this  paper,  we  provide  an 
approximation  to  /*(X«  A)  which  is  both  easy  to  compute  and  accurate  over  a  wide 
range  of  parameter  values.  In  typical  applications,  A  is  an  orthant  {x:  x  ^  c; 

i«1,...,n),  a  union  of  orthants  {x:  j£l(xk  £  c)  «  p},  a  cube  (x:  -c  $  x.  ^  c;  i«1 n},  or 

2 


a  sphere  (x:  jExfc2  £  c). 


The  problem  of  evaluating  P(X(  A)  arises  in  many  contexts.  In  genetic  models  of 
the  transmission  of  traits,  this  probability  is  related  to  an  individual's  risk  of 
contracting  a  disease:  see  Rice,  et  at. [24],  Henery[8],  and  Curnow[2].  In 
geostatistics,  it  is  related  to  the  probab:lity  of  correct  classification  in  indicator 
kriging:  see  Journel[14],  and  Solow[28].  It  is  useful  in  the  construction  of 
simultaneous  confioence  intervals:  see  MillerC  19],  and  Uusipaikka[32].  It  also 
appears  in  the  middle  of  some  time  series  or  regression  problems:  see  especially 
Keenan[  15],  who  shows  that  orthant  probabilities  are  prediction  error  probabilities 
and  transition  probabilities  for  certain  binary  time  series,  and  Keener  and 
Waldman!  16],  who  need  to  evaluate  orthant  probabilities  (for  dimensions  as  large  as 
35)  for  computing  the  likelihood  function  for  rank-censored  data.  Our  restriction  on 
the  set  A  may  seem  quite  severe;  however,  this  formulation  is  already  useful  in  the 
applications  just  mentioned.  In  fact,  the  probability  content  of  the  positive  orthant 
has  been  the  subject  of  many  papers:  see  Martynov!  17],  and  Moran[21].  Also,  the 
probability  content  of  the  sphere  gives  the  distribution  of  a  weighted  sum  of 
independent  X2  variates:  see  Solomon  and  Stephens[27],  and  Imhof[l0],  for  other 
approaches. 


If  a  multivariate  normal  probability  were  required  only  infrequently,  then  ordinary 
simulation  or  numerical  quadrature  would  often  be  adequate.  However,  in  many 
applications  (see  especially  [15]  and  [16])  the  computation  of  such  probabilities  is 
just  a  small  component  of  a  much  larger  problem.  This  larger  problem  typically  — 
involves  the  estimation  of  parameters,  and  uses  an  iterative  algorithm;  thus  a-  ~ 
subroutine  which  evaluates  multivariate  probabilities  is  repeatedly  called,  and  so  □ 

speed  and  accuracy  are  important  considerations. 


ru 


This  problem  is  an  old  one,  and  various  approaches  have  been  proposed;  we  now 


briefly  review  them.  Numerical  integration  and  simulation  in  high  dimensions  are  'L. 


□ 


known  to  be  quite  slow.  Recently,  however,  Schervish[26]  wrote  a  program  for  'y  Col®s_ 

tnd/o  r 


■r 
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computing  rectangular  probabilities,  and  tested  it  for  dimensions  n^B,  end  s  p. 
Simulation  is,  of  course,  always  available  in  its  simplest  form,  but  standard  variance 
reduction  techniques  are  hard  to  obtain;  a  notable  exception  is  the  recent  work  of 
Moran[2l],  in  which  he  gave  a  clever  method  for  estimating  the  probability  of  the 
positive  orthant,  and  provided  a  control  variate  for  the  estimator.  Plackett[23]  gave 
a  reduction  formula,  but  it  seems  practical  only  for  r>^5.  The  study  of  probability 
inequalities  (see  Tong[31],  and  Eaton[5])  has  provided  many  interesting  insights  and 
useful  results,  but  the  inequalities  themselves  often  yield  poor  approximations;  see 
Section  5  below  for  further  discussion.  The  evaluation  of  multivariate  normal 
probabilities  is  closely  related  to  the  determination  of  volumes  of  sets  on  the 
surface  of  the  unit  ball  in  ]ftn,  but  this  relationship  has  only  occasionally  been  used: 
see  for  instance  the  work  of  Abrahamson[  1  ],  and  Ruben[25];  there  is  a  recent 
revival  of  geometrical  methods  in  various  contexts,  as  seen  in  the  work  of  Diaconis 
and  Efron[4],  Johansen  and  Johnstone!  12],  and  Naiman[22].  Next,  expansions  such 
as  the  tetrachoric  series  have  been  tried  in  several  cases,  but  difficulties  here  include 
their  slow  convergence  or  even  divergence  over  much  of  the  parameter  space:  see 
Harris  and  Soms[7],  and  Moran[20].  Finally,  sometimes  a  direct  approach  is 
rewarding.  For  instance,  if  p.^  =  p  £  0,  or  if  />..*«.«.  where  -1<a.<1,  then  many 
relevant  computations  reduce  to  the  evaluation  of  one-dimensional  integrals,  which 
are  easy  to  do:  see  Steck[29],  and  Curnow  and  Dunnett[3]. 

In  this  paper,  we  will  exploit  the  symmetry  of  A  to  provide  an  approximation  to 
/*(X »A)  that  is  expressible  in  terms  of  one-dimensional  integrals;  such  integrals  are 
easily  evaluated  by  computer.  Some  of  the  methods  mentioned  above  come  into 
play:  for  example,  we  will  use  the  case  /> ..  =  p,  about  which  we  will  expand  a  Taylor 
series  which  turns  out  to  be  an  (integrated)  Gram-Charlier  series,  and  we  will  use 
Schervish's  program  for  comparison.  After  establishing  our  notation  in  Section  2,  we 
introduce  a  simple  approximation  in  Section  3.  We  compute  and  study  correction 
terms  in  Section  4,  do  two  examples  in  some  detail  in  Section  5,  and  summarize  our 
numerical  work  there.  We  then  conclude  with  a  brief  discussion  in  Section  6. 

2.  NOTATION  AND  PRELIMINARIES 

Let  X  *  (X,,...,X  )'  be  a  standardized  nonsingular  multivariate  normal  random  vector 

1  r> 

whose  distribution  is  denoted  N  (0,R);  thus  E(X)  «  0,  and  var(X)  *  R  «  (/>  )  with  p*  1 

n  r  ij  r  u 

for  all  i.  The  density  of  X  is  denoted  f  (x;R);  the  univariate  normal  distribution  is 
denoted  $00  «  ^(X^x)  and  its  density  is  fix).  Let  d  *  n(n-1)/2,  and  string  out  the 
correlations  into  a  vector  /*“(^12'/,1 IRd;  B,s0'  ,et  ~p  be  the  8ver89e 
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of  the  components  of  />,  end  let  p  «  { p t  ]Rd.  Ea  is  an  n-by-n  equicorrelation 
matrix  with  parameter  «,  where  -<n-1)“ 1  <  «  <  1.  Pn  is  the  set  of  n-by-n  permutation 
matrices,  an  element  of  which  is  denoted  >;  II  Is  a  random  matrix  which  is 
independent  of  X  and  is  uniformly  distributed  over  P  .  The  permutation-symmetric 
set  A<IRn  satisfies  A  *  »A  «  {»a:  a<A}  for  all  v  t  P  .  Our  main  example  below 
involves  the  positive  orthant  Q  ■  {x:  x,  ^  0;  k*1....,n)  *  {x;  x  ^  0}.  Unless  stated 
otherwise,  all  summations  are  over  the  entire  range  of  the  index  variable:  thus,  for 
instance,  I  is  a  summation  over  all  of  P  .  Finally,  let 

"  n 

(1)  h  (p)  =  P(X«A). 

We  suppress  the  dependence  of  h^  on  A. 

To  derive  his  reduction  formula,  Plackett[233  proved  and  used  the  following  identity 
for  the  multivariate  normal: 

(2)  fc/a/»  >  ^n(x;R)  =  <a2/ax.ax.)  ^(x;R). 

This  identity  has  been  used  to  establish  probability  inequalities  and  monotonicity 
results  for  multivariate  normal  probabilities  (see  Tong[31]).  We  use  it  repeatedly 
below:  it  simplifies  many  computations,  and  is  a  natural  way  of  deriving  the 
multivariate  Gram-Charlier  series  in  Section  4. 

The  equicorrelated  case  is  important  in  our  discussion  below,  so  we  describe  it 

now.  If  Z  =  <Z„„..,Z  )'  is  a  N  (0,1)  vector,  Z  is  a  standard  normal  independent  of  Z, 
inn  o 

e  =  (1,...,1)’  t  lRn,  and  «>0,  then  V  *  (1-«)1/2Z  ♦  «1/2Z  e  is  a  N  (0,E  )  variate.  Upon 
conditioning  on  Z  ,  we  get  the  single  integral 

rA  -  t«1/2e 

P(Z  t  ( - m)  )  f[t)  dt. 

1 17 

-00  (1  -«)17^ 

When  A  is  an  orthant  or  a  cube,  the  integrand  in  (3)  is  just  a  product  of  one- 
dimensional  normal  marginal  probabilities,  and  if  A  is  a  sphere,  it  is  a  noncentral  A 
probability.  Thus,  for  many  cases  of  interest,  the  right  side  of  (3)  is  easily 
evaluated.  The  analogous  formula  for  «< 0  is  more  involved,  but  still  tractable:  see 
Steck[29].  Also,  moments  of  the  form 

(4)  EIV,V  .  .  .  V  *n  KViAlJ 

1  n 


are  needed  below,  and  they  can  be  expressed  as  a  single  integral  just  like  (3).  The 
same  argument  applies  to  the  case  />..»«.«.,  with  -1  <  «.  <  1  for  all  i,  which  is 
generated  by  V  «  (1-a  2)1/2Z  ♦  a  2  ,  i«1,...,n. 

3.  A  SYMMETRY  ARGUMENT 

Since  A  *  »A  for  all  kP,  P{7Lt A)  =  P(srX« A)  for  all  P  ,  and  P{7L* A)  *  PfriXrA) 

n  n 

also.  Consider  I1X;  it  has  exchangeable  components  and  it  is  a  scale  mixture  of 
normals  with  density 


(5)  f  (x;R)  *  (nl)_12:  f  (x;»Rff'). 

n  —  »  n  - 

Its  first  two  moments  are  E(I1X)  *  0  and  var(IIX)  *  (n!)~ 1  ^  ir R*r '  *  E-.  _  We  get  our 
simplest  approximation  by  fitting  to  fn  the  normal  density  ^(x;E^),  which  shares  the 
same  first  two  moments:  thcs. 


(6)  h  (/>)  =  h  (p). 

This  approximation  has  several  appealing  properties.  First,  it  is  easily  evaluated  for 

many  cases  by  the  argument  leading  to  (3).  Second,  there  is  the  following  heuristic 

argument.  Let  P.  be  the  set  of  d-by-d  permutation  matrices,  and  let  P  '  be  the 

subset  of  it  induced  by  the  correlation  vectors  of  (»X:  arP^}.  The  points  {»/>: 

vfP  '}  are  the  vertices  of  a  regular  polygon  centered  at  p.  Since  h  has  the  same 

value  at  each  vertex  of  this  polygon,  and  since  it  is  a  smooth  function  of  p.  its 

value  at  the  center  of  the  polygon  should  not  be  far  from  its  value  at  a  vertex. 

Third,  it  comes  from  the  least  squares  fit  of  the  equicorrelation  matrices  to  R;  that 

is,  it  minimizes  Z(/>  -a)2  over  a.  And  fourth,  a  Taylor  expansion  of  h  about 
i.j  ij  n 

«=(<*,...,«)'  gives 


(7)  h  (p)  -  h  (c)  +  Dh  (a)'(p-a)  *  remainder, 

n  .  n  -  n  ~  •  - 


where  D  denotes  the  gradient.  Each  component  of  Dh  <«)  is,  by  Plackett's  identity 

n  — 
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“  3a  (92/3x|3xj)  d!  "  3a  02/3x13x2>  *„<£*«*  d5* 

The  interchange  of  the  order  of  differentiation  and  integration  in  (8)  is  easily 
justified,  and  the  last  equation  there  comes  from  the  symmetry  of  A.  The  linear  term 
in  (7)  is  thus  zero  if  o*p,  so  that  (6)  is  also  a  good  first  order  approximation. 

Approximation  <6)  <s  based  on  only  two  moments,  so  it  is  not  expected  to  be 
accurate  for  extreme  tail  probabilities,  for  example  for  />(X.  2  c;  all  i)  for  large  c; 
see  Steck[30]  and  lyengar[11]  for  a  treatment  of  such  cases.  Also,  it  is  clear  from 
the  heuristic  argument  above  that  the  approximation  should  improve  as  Z(p.-p)2 
decreases.  Thus,  we  need  higher  moments  of  FIX  to  provide  correction  terms,  sc  we 
now  turn  to  the  higher  order  terms  of  the  Taylor  expansion  (7). 

4.  CORRECTION  TERMS 


To  write  the  Taylor  expansion  of  h  (/>)  about  h  (p).  we  need  additional  notation:  let 

n  -  n  ~ 

k  =  (k)2 . kn_ln)  consist  of  nonnegative  integers  k.  .  for  i<j,  and  let  |k!  ■  k1  2  ♦  ...  ♦ 

*  2 !  k  ! .  Next, 


n-1.n' 

let 


;  let  m  =  (m,  ,...,m  ).  where  m.  =  Z  k  «•  Z 

1  n  1  .-.1,1  ... 

i<i  .>j 


j.' 


and  note  that 


!  m ! 


(91  /  *  />,  J1''2-  '  -Vm"-''* 

(10)  Dkhn(^)  *  0P/3/»12k  1.2.  .  ^/>n_1nk n-l.n)  h^).  for  !k|  *  p; 

(11)  D">n(x;R)  =  O^/dx^i.  .  ^xnmn)  ^n(x;R),  for  imj  =  2p. 


(12) 


C(p.k)  =  p!  /  (k,  J.  .  .k  ,  !),  for  Ik  I 

«  1,4  n- 1  ,n  • 


p. 


We  then  have 


(13) 


h„'f>  ■ 


Z  C(p,k)  (p  -  ~p)k  Dkh  (/>)  +  remainder. 

Ikl-p  ‘  “ 


Consider  the  integrands  in  (13):  on  the  left  side,  we  have  jMx;R);  by  changing  the 
order  of  differentiation  and  integration  on  the  right  side  and  by  repeatedly  using 
Plackett's  identity  (2)  there,  we  get 


(14) 


t&rt  *  2v$Pfr1 


Z 

Ik  I  *p 


C(p.k)  (p  -  p)k 


*  remainder 
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«  I^Jpir 1  2  ^  C(p,k)  (/>  -  p)k  H^txjE^)  ♦  remainder. 

where  HMxjE^  is  the  m111  order  Hermite  polynomial.  Thus,  the  Taylor  expansion  (13) 
yields  the  Gram-Charlier  series  expansion  of  jMx;R)  about  ^n(x;E ^);  see  Johnson  and 
Kotz[13]  and  Mihaila[18]  for  a  discussion  of  the  multivariate  Gram-Charlier  series 
and  Hermite  polynomials.  Of  course,  we  could  also  interpret  the  integrands  of  (13) 
as  the  Gram-Charlier  expansion  of  the  mixture  of  normals  ^n(x;R)  about  ^(x;E^;  we 
get  a  different  expansion,  but  either  interpretation  gives  the  same  correction  terms 
for  the  approximation  (6). 


For  our  purposes,  the  important  feature  of  (13)  is  that  each  Dkhn(^)  can  be  expressed 
as  a  one-dimensional  integral;  this  is  because  of  the  Hermite  polynomials  of  (14)  and 
the  argument  leading  to  (3)  and  (4).  It  is  easy  to  see  that  the  p**1  term  of  (13) 
requires  (on  the  order  of)  pd  additions,  where  d  =  n(n-1)/2.  Note  that  no  matrix 
inversion  is  necessary  to  evaluate  the  correction  terms.  In  the  next  section,  we  give 
details  of  one  special  case  to  illustrate  our  use  of  (13). 


The  convergence  of  (13)  is  a  difficult  issue  in  general.  One  special  case  that  is 
tractable  is  the  trivariate  positive  orthant,  Q  ,  for  which  h  (/»)  *  1/2  -  (1/4») 
<cos_1(/>12)  *  cos_1{/>13)  ♦  cos ~}[p23)).  Here,  (13)  comes  from  the  expansion  of  cos-1 
about  ~p.  For  any  fixed  />  £  0  or  ~p  £  3-2(3 1  /2>,  (13)  is  convergent  for  all  choices  of 
p  leading  to  that  />  and  yielding  a  positive  definite  matrix  R.  However,  when  ~p  t 
(3-2(3 1/2),0),  (13)  is  not  always  convergent:  take  for  instance  p  *  (0.45,-0.60.-0.75)';  the 
radius  of  convergence  for  cos-1  about  ~p  *  -0.30  is  0.70,  and  />12  *  0.45  lies  beyond 
that.  We  will  rigorously  study  the  convergence  of  (13)  and  (14)  elsewhere;  here,  we 
depend  upon  numerical  examples  to  assess  the  performance  our  approximation. 


5.  EXAMPLES 


Our  first  example  deals  with  the  positive  orthant  Q  =  (x:  x  >  0).  We  give  explicit 
expressions  for  the  first  two  terms  of  the  Taylor  expansion  (13),  and  describe  their 
numerical  evaluation.  The  algebra  here  is  straightforward  but  rather  tedious,  so  we 
suggest  the  use  of  MACSYMA  for  other  applications.  We  only  do  the  case  0  here. 

% 

The  first  term  of  the  Taylor  expansion  is 

(15)  hn(^)  «  ♦(tr1/2)n  ^(t)  dt, 

n "  -oo 
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where  r  >  pl{\-p).  This  expression  cen  be  accuretefy  evaluated  using  a  Gaussian 
quadrature  formula  and  Hill’s  [9]  algorithm  for  computing  ♦.  For  n^S,  quadrature  is 
not  necessary;  this  is  because  (15)  reduces  to  1/2,  1/2  -  (1/2»)cos ~'(p),  and  1/2  - 
(3l4w)co*~'(p)  for  n  *  1,  2,  and  3,  respectively.  For  n*4,  a  simple  formula  is  not 
available,  but  Steck[29]  provides  an  accurate  approximation;  and  for  n=5.  we  have 

(16)  h  (p)  *  — ♦  -Lcos_1(^)  ♦  Lh<* 

We  need  (15)  for  n^O  below,  so  define  h  _^)  s  0  for  n=1,2,...,  and  h0<^)  *  1.  We  also 
need  the  following  restricted  moments: 

(17)  f  (e;a . a)  =  E[V,*i  .  .  .V*t  I(V<Q  )]. 

nit  1  t  n 


where  V  is  a  Nn(0.Ea)  variate,  a  i  .  .  .  ^a^l,  and  l^t^n.  Using  the  argument  leading 
to  (3),  we  can  rewrite  the  integral  in  (17)  as  a  one-dimensional  integral:  for  instance, 

(18)  fnU;1)  *  (1  -«)1/2  ♦  ti*(t«)]  *(t*)"'1  fit)  dt, 

?l  -oo 

where  o  *  («/ti-«))  .  For  smail  values  of  n,  this  integral  can  be  evaluated 

explicitly;  else,  numerical  quadrature  is  required. 

Next,  the  quadratic  form  in  (13)  has  an  expansion  whose  terms  are  proportional  to 

(19)  (p  -p)(p^-p)  JQ  04/ax.axj^xkax|)  ?n(x;E^)  dx, 

with  i<j  and  k<l.  Such  a  set  {i.j,k,l}  has  two,  three,  or  four  distinct  elements.  If 
there  are  four,  then  the  integral  in  (19)  is 


(20) 


o;e>) 


dx. 


where  a  *  p/^^*4p);  the  integral  here  is  evaluated  in  the  same  way  as  (15).  When 
(i.j.k.l}  has  three  distinct  elements,  the  integral  in  (19)  is 


(21)  (n-3)  f3( 0;E^)  A (p)  f„_3<Al). 


where 
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(22)  A •  ~p  iU3p)'12  I  (1-(n-1)>)  (V\-/MU2~p))'12, 

and  fi*pl(^3fi.  For  n^6,  the  expectation  f  in  (21)  can  be  evaluated  explicitly  in 

n- j 

terms  of  elementary  functions;  we  omit  those  details.  And  when  {i,j,k,l}  has  two 
distinct  elements,  the  integral  in  (19)  is 

(23)  *2<0;E^>  [B(p)  (fn_2(y;2Hn-3)fn_2(r;1.D)  ♦  C(/JQ  *n_2<x;Ey)  dx], 
where 

(24)  BCp)  *  (n-2)<1*2/>)  /  (1-/>2)  (T(n-1)^>2, 

(25)  C{p)  *  ~p  I  {1-~p)  (1*(n-1)>), 

and  y-~pH“\*2~p);  the  expectations  f^_2  here  can  be  explicitly  evaluated  without 
numerical  quadrature  if  rv^5.  Thus,  we  can  use  the  expressions  (20),(21),  and  (23)  to 
get  the  first  correction  term.  We  omit  the  next  correction  term  which  involves  the 
third  derivatives  of  h  (7>). 

n  C 

We  now  turn  to  our  numerical  work.  For  dimensions  n  *  3,  4,  and  5,  we  generated 
correlation  matrices  R  with  p  ^  0  for  all  i,j.  For  each  case,  we  computed  the  true 
probability  of  the  positive  orthant  using  the  program  of  Schervish.  We  also 
computed  approximation  (6)  and  corrections  from  the  second  and  third  derivatives  of 
(13).  We  then  plotted  the  relative  error  of  the  three  approximations  against  the 
"variance  in  p":  (n)_1  X  (p  -  ~p)2\  these  plots  are  shown  in  Figs.  1-9. 

For  n=3,  the  relative  error  decreases  considerably  as  we  add  higher  order  terms  of 
(13).  Two  parts  of  the  scatterplot  are  clearly  separated  in  Figs.  1  and  3  (and 
somewhat  less  clearly  in  Fig. 2):  the  upper  part  corresponds  to  "extreme"  p  of  the 
form  (0.0,/>13.0.9),  and  the  lower  part  corresponds  to  less  extreme  p.  For  n  *  4,  the 
relative  error  of  the  third-order  approximation  is  often  smaller  than  that  of  the 
second-order  one,  although  the  difference  between  the  two  is  not  great;  thus,  the 
range  of  relative  errors  shown  in  Figs.  5  and  6  are  the  same.  For  n  *  5,  this 
improvement  is  somewhat  greater.  A  conclusion  from  our  work  is  that  (for 
correlation  matrices  of  this  kind)  three  terms  of  the  series  (13)  seem  to  be  enough  to 
give  a  relative  error  of  less  than  about  seven  percent. 
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Next,  we  turn  to  the  speed  of  our  algorithm.  In  Table  1,  we  compare  the  time 
required  (Tf)  to  compute  the  first  three  terms  of  our  approximating  series  with  the 
lime  required  (T^)  to  compute  the  exact  probability  (demanding  three-place  accuracy) 
using  Schervish's  program.  Alt  times  are  given  in  seconds;  they  are  the  average  of 
five  runs  on  a  VAX  750.  For  all  cases  cited  in  Table  1,  the  relative  error  of  our 
approximation  was  less  than  0.03.  Notice  that  Schervish's  program  is  much  faster 
when  p  is  close  to  />  than  when  it  is  not.  The  reduction  in  the  time  needed  is 
especially  evident  for  n  *  5  when  the  correlation  matrix  is  not  too  close  to  the 
equicorrelated  case. 

Our  second  example  involves  exceedance  probabilities.  For  a  fixed  constant  c^O, 
let  S  «  ZI(X  ^  c),  and  let  p  ,(/>)  =  P(S  =k).  Then 

n  x  k  rn,k  -  n 

(26)  pn<£)  =  (P  nJp) . pjp)) 


is  the  exceedance  distribution  which  we  approximate  by  p Jp).  Note  that  the  mean  of 
the  exceedance  distribution  and  that  of  its  approximation  are  the  same:  E(S  :p)  - 

n  — 

E(S  ;p)  *  ni’(-c).  We  now  study  the  second  moment  of  the  exceedance  distribution 

n  — 

and  that  of  its  approximation;  this  study  will  yield  new  results  about  our 
approximation  of  orthant  probabilities.  The  variance  is  given  by  var(S  ;/>)  =  n$(-c)  - 

n  — 

n^-c)2  ♦  X  m  £  c,  X  £  c),  which  depends  only  on  bivariate  quadrant  probabilities. 
We  need  the  following  lemma,  the  proof  of  which  is  an  easy  application  of 
Plackett's  identity. 

LEMMA  1:  Let  V  be  a  N2(0.R)  variate  with  p  J2  *  r.  Then  F(r)  *  P(M  ^  ^  c,  V2  ^  c)  is 

convex  for  r  in  (0,1);  for  c  £  21/2  -  1.  F  is  convex  in  (-1,1);  and  for  c=0,  F  is  concave 

in  (-1,0). 

This  lemma  yields  the  following 

PROPOSITION  1:  If  (i)  p ^  £  0  for  all  i  and  j,  or  (ii)  if  c  £  21/2  -  1,  then  var(S ^p)  £ 

var(S n;/>);  if  (iii)  c=0  and  p  £  0  for  all  i  and  j,  then  var(Sn;/>)  £  var(Sn;^). 

Proof:  The  function  g (/>)  «  var(S  ;p)  is  a  symmetric  function  of  the  components  of 

—  n  - 

p.  Under  (i)  or  (ii),  g  is  Schur  convex  function  of  p  (see  Tong[31,p.106]);  and  under 
(iii),  g  is  a  Schur  concave  function  of  p.  Since  p  majorizes  p,  the  result  follows. 

Now  consider  the  case  of  the  positive  orthant  (c*0,  and  Sn=n)  and  suppose  that  p 
^  0  for  all  i  and  j;  let  «  *  (« . a),  where  a  «  min  { p  )  and  /  «  where  fl  * 
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max  From  the  variance  inequality  in  Proposition  1.  it  is  tempting  to  conjecture 

that 

(27)  p  (p)  £  p  Ip). 

n.n  -  an  - 

This  is  indeed  true  when  n=3  (because  -cos'1  is  convex  in  (0,1)),  and  it  represents  an 
improvement  over  the  bounds  derived  from  the  result  of  Slepian:  p  (a)  ^  p  (p)  < 

an  —  an  tm 

p  (fi)  (S8e  Tong[3l,p.10  and  p.169]).  Examples  show  that  inequality  (27)  is  not  true 

an  — 

for  larger  n.  However,  Figures  10-12  show  that  the  right  side  of  (27)  is  a  much  better 
approximation  to  the  positive  orthant  probability  than  are  the  Slepian  bounds:  for  the 
same  correlation  matrices  studied  above,  these  Figures  compare  the  relative  error  of 
the  Slepian  bounds  with  the  relative  error  of  our  approximation  (6);  our  approximation 
was  better  in  every  case  (except,  of  course,  when  the  the  matrix  was  already 
equicorrelated,  in  which  case  all  three  coincided  with  the  true  probability.) 

The  same  method  gives  similar  results  for  the  cube:  T  *  ?  I(|Xt|  <  c),  q  ,  s 

n  pc.  k  n,k 

/^T^k).  The  analogs  of  Lemma  1  and  Proposition  1  come  from  the  behavior  of  G(r) 
=  P[\My  \  s  c.:V2!  <;  c). 

6.  DISCUSSION 

The  problem  of  evaluating  multivariate  normal  probabilities  is  a  difficult  one,  and  it 
is  likely  that  there  is  no  panacea;  that  is,  an  approximation  which  is  tailored  to  work 

well  for  one  set  of  parameter  values  will  probably  be  inadequate  for  another.  Thus, 

many  methods  have  been  proposed  in  the  literature.  Our  contribution  to  this 
literature  is  to  show  (from  our  theoretical  and  numerical  work  above)  that  (13) 

provides  a  good  approximation  to  probabilities  of  permutation-symmetric  regions,  and 
that  it  is  easily  evaluated.  One  drawback  of  our  proposal  is  that  useful  error  bounds 
for  our  approximation  are  not  yet  available;  here,  we  depend  upon  numerical  work  to 
assess  the  error.  Our  methods  can  also  be  applied  to  the  problem  of  evaluation 
Ef(X),  where  f(x)  =  f(irx)  for  all  v « P  ,  and  also  to  random  variables  with  other 

n 

elliptically  contoured  distributions. 
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Figure  11. 
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Figure  12. 
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Table  1 
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1 
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( »2» *2  f • 2  9 *2# • 2  f *2) 

0.81 

1.51 
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(•2#  »4| #4» • 6 > • 8 1 *8) 
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14.60 
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(.0,.0f.0,.2,.2,.2,.4,.4,.4,.4) 
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